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Local phonon density of states in an elastic substrate 
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The local, eigenfunction-weighted acoustic phonon density of states (DOS) tensor is calculated for 
a model substrate consisting of a semi-infinite isotropic elastic continuum with a stress-free surface. 
On the surface, the local DOS is proportional to the square of the frequency uj, as for the three- 
dimensional Debye model, but with a constant of proportionality that is considerably enhanced 
compared to the Debye value, a consequence of the Rayleigh surface modes. The local DOS tensor 
at the surface is also anisotropic, as expected. Inside the substrate the local DOS is both spatially 
anisotropic and non-quadratic in frequency. However, at depths large compared with vi/u, where 
vi is the bulk longitudinal sound velocity, the local DOS approaches the isotropic Debye value. The 
results are applied to a Si substrate. 

PACS numbers: 68.35.Ja, 05.70.Np, 46.70.-p 
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INTRODUCTION AND ELASTIC 
SUBSTRATE MODEL 



A variety of surface and near-surface thermodynamic 
and lattice dynamical phenomena are governed by the 
phonon density of states (DOS), the number of vibra- 
tional modes per unit frequency, per unit volume of solid. 
In an inhomogeneous system, however, it is possible to 
distinguish between the "global" DOS, which makes ref- 
erence to only to the vibrational-mode frequencies of the 
entire system, and a local, eigenfunction-weighted DOS. 
In systems with a small mechanical surface-to- volume ra- 
tio, the volume including all parts of the system and envi- 
ronment participating in the vibrations, the global DOS 
is essentially unchanged from that of a bulk system. The 
eigenfunction-weighted DOS is interesting both because 
it is relevant to a variety of physical phenomena, and 
because it contains local, position-dependent mechanical 
and vibrational information about the system. 

In this paper, we calculate the local phonon DOS ten- 
sor for a model substrate consisting of a semi-infinite 
isotropic elastic continuum with a stress-free surface, 
appropriate for low-frequency acoustic phonons in real 
solids with an exposed clean surface. The elastic contin- 
uum is also assumed to be linear and nonpolar. The free 
surface is taken to lie in the xy plane, with the substrate 
occupying the z > half-space. Periodic boundary con- 
ditions are applied in the lateral directions, on an area 
A. Normalized vibrational eigenfunctions have been ob- 
tained for this geometry in a classic paper by EzawaJ- 



where u(r, t) is the displacement field, and where u^ = 
(diUj + djUi)/2 is the strain tensor, p is the mass density, 
and A and /i are the standard Lame coefficients, related 
to the bulk longitudinal and transverse sound velocities 
vi and vt according to vi — \/{X + 2ii)/p and vt = \/fi/p- 
The equation of motion following from Eq. Q is 



dfu^ 



vfViV ■ u) 



Uj V X V X u 



(2) 



It will also be useful to define an elastic stress tensor 
Tij from the continuity equation 9(11^ + djTij — for 
momentum density 11 = pdtu. In an isotropic elastic 
medium, it follows from Eq. (O that 



T ■ — 



-A(V • u)(5y — 2p,Ui 



-Cijkl Ukl, 



where 



Ctjki = XS^jSki + p.{5ik5ji + SuSjk) 



(3) 



(4) 



is the corresponding elastic tensor. The boundary condi- 
tion at the z — Q exposed planar surface S is 



T,,(r) = 0, reS, 



x,y,z. 



(5) 



To quantize the elastic waves we require the displace- 
ment field and momentum density to satisfy the canonical 
commutation relation 



K(r),n,(r')]-iM,,<5(r-r'). 



(6) 



Therefore, we can expand the displacement field in a ba- 
sis of bosonic creation and annihilation operators accord- 
ing to 



II. LOCAL PHONON DOS TENSOR 

We begin by defining the local DOS tensor gij{r,uj), 
the quantity to be calculated in this paper. The La- 
grangian density for a linear, nonpolar, isotropic elastic 
continuum is 



C = hp{dtu) 



jAujj 



f^ufj, 



"(r) = E 



2puj„ 



a„f„(r)-|-aj,f;;(r) 



(7) 



where the fn(r) are vibrational eigenfunctions, defined 
to be time-periodic solutions of the elasticity equation 
(0), in the presence of the stated boundary conditions, 
and normalized according to J(Pr {* ■ {„' = <5„„'. Here 
(1) the integration is over the volume V of the system. The 



vibrational eigenfunctions are assumed to be complete, 
sothatS„/^(r)/r(r')=<5,,<5(r-r').^ 

The quantity we shall calculate in this paper ia^ 

g,,{r,cu) ^ ^ 5(a; - c^„) /^(r) [fiir)]*, (8) 

n 

where i and j are Cartesian tensor indices. This quantity 
weights each vibrational mode n by the square of the 
amplitude of the eigenfunction at position r. The trace of 
gij{r,Lu) also contains physically interesting information, 
namely 



where m = (1, ti, t2) labels the longitudinal and two trans- 
verse branches, the sound velocities are the bulk values 
quoted in Sec. E] and the emk are unit polarization vec- 
tors forming an orthogonal basis, eik is orietnted in the 
k direction, and et^k and etjk are perpendicular to eik 
and to each other. 

In this case, Eq. ^ yields 



5y(r,^) 



-y 

mk 



b{{jJ - UJ,n\^ i 



mk 



(13) 



Tr g,j{i 






9ii{^,^) 



>(^-^„)|f„(r)|2 (9) 



Now, averaging the trace of 5y(r, lu) over the system vol- 
ume V leads to 



1 /" 1 .. 



(10) 



which is the ordinary intensive global DOS, i.e., the num- 
ber of vibrational modes per unit frequency per volume. 
In an isotropic system with translational invariance, 
the local DOS is independent of position, and Eq. (0) 
simplifies to 



gij{r,uj) = ^^(5(cj-w„), 

n 



(11) 



where a three-dimensional system is assumed. Because 
there is no position dependence here, the trace of this 
quantity yields the global DOS [the right-hand-side of 
Eq. ^]. 

Is the local DOS defined Eq. © a quantum mechanical 
quantity? One the one hand, it does not involve h, and 
in that sense it is not quantum mechanical. On the other 
hand, however, the normalization assumed for the f„(r) 
is of quantum-mechanical origin: Assuming the conven- 
tional bosonic commutation relations for the a„ and ajj 
appearing in the expansion of Eq. |(7J), the normalization 
condition on the vibrational eigenfunctions is required for 
Eq. © to be satisfied.^ Therefore, the local DOS defined 
Eq. (jHl naturally shows up in quantum mechanical cal- 
culations of the vibrational properties of inhomogeneous 
systems. 



III. LOCAL DOS IN THE 
THREE-DIMENSIONAL BULK LIMIT 

In what follows it will be useful to recall the well-known 
Debye DOS formula. The spectrum and normalized vi- 
brational eigenfunctions for a three-dimensional isotropic 
elastic continuum, with periodic boundary conditions ap- 
plied to the surface of a three-dimensional volume V, are 



^mk = Vr, 



and fmk(r) 



o^kr 



W 



^TTlki 



(12) 



The fact that Wmk is independent of the direction of k 
allows the product of polarization vectors to be averaged 
over the unit sphere, leading to Sij/D in D dimensions. 
We then obtain 



?TT,k 



(14) 



an example of the form given in Eq. Hll|) . Therefore, in 
the V — > oo limit. 



Sij Lu'W I 2 



ff.,(r,c.) = — X— (-3+-3 



(15) 



Eq. (|15|l is the conventional three-dimensional Debye 
DOS tensor for acoustic phonons. Of course, in a real 
solid the continuum approximation breaks down at short 
wavelengths and the spectrum is cut off at high frequen- 
cies. 



IV. LOCAL DOS IN SEMI-INFINITE ELASTIC 
SUBSTRATE WITH A STRESS-FREE SURFACE 



The evaluation of the local DOS for our substrate 
model requires the vibrational eigenfunctions of a semi- 
infinite isotropic elastic continuum with a stress-free sur- 
face, which have been obtained by Ezawai The modes 
are labeled by a branch index m, taking the five values 
SH, -|-, — , 0, and R, by a two-dimensional wavevector K 
in the plane defined by the surface, and by a parameter 
c with the dimensions of velocity that is continuous for 
all branches except the Rayleigh branch m = H. The 
range of the parameter c depends on the branch m, and 
is given in Tabled The frequency of mode (to,K,c) is 
^mKc — cK, where K = |K| is the magnitude of the 
two-dimensional wavevector. As stated above, we shall 
assume periodic boundary conditions in the x and y di- 
rection applied to a finite area AS 



m 


range of c 


SH 


[Vt , CXj] 


± 


[1)1,00] 





[«t,ul] 


R 


cr (discrete) 



TABLE I: Values of the parameter c for the five branches of 
vibrational modes of a semi-infinite substrate. The Rayleigh 
branch m — K has a single value of c, which we denote by cr. 



The definition of the local DOS given in Eq. applies to a system with a discrete spectrum. The generalization 
of Eq. (jSJ to our model, which has a mixed discrete-continuous spectrum, is 



(r, Lo)=Y^ S{u cnK) fi,^{r) [f^^ir)]* + J] J] / dc 5(c. - cK) fUdr) ifUdr^)] 

K K m^Ry^" 



(16) 



The first term accounts for the Rayleigh surface wave branch to = R, and these vibrational mode eigenfunctions are 
normalized as in Sec.^ The value of c for the Rayleigh branch, cr, is calculated below. The second term in Eq. H16(l 
accounts for the four remaining branches m = SH, -I-,— , and 0, having continuous values of c. Their eigenfunctions 
are instead normalized according to JdPr f^j^^ ■ fm'K'c' = ^mm' (^kk' S{c -— c'). Furthermore, Tm is the integration 
domain given in Tabled 



We turn now to an evaluation of Eq. (|16|l . The con- 
tribution from each branch in Table Q] will be discussed 
separately in the sections to follow. 



A. SH branch 

First we consider the m — SH branch, for whichi 



ISH 



/ 2cK 



cos{fiKz) Gz X gk e 



iK r 



where 



I3{c) ^ Vic/vtr - 1, 



(18) 



and where ex = Ji./ K is a unit vector in the K direction. 
This mode is polarized in the xy plane, perpendicular to 
K. The contribution made by this branch to the local 
DOS is 



Hi] 



(r,^) 



2uj 



2a Z^ 

t^ K -^^t 



dc 



cos^ [to [iz/c) 



X 5{ijj-cK) {e^xe^y{e^-Ke]:^y . (19) 



Taking the IK ^ 00 limit, and making use of the two- 
dimensional DOS formula 






K 

leads to 



4^ 



S{lo - cK) 



27rc2' 



(20) 



JSH) 



(r,c.) 



T^^VI j^^ 



dc 



cos^ {ijjj3z/c) 



(ezXeK)''(e2XeK)-'', 
(21) 



where 



//70 
2:^ (ezXeK)'(ezXeK)^' (22) 

denotes a two-dimensional average over the directions of 
K. Finally, using 



5iA 



(ezXeK)*(ezXeKy = ^-(^ 



(17) leads to 



9ij (r, uj) 



^^^5..(l-<5..)/SH(f) 



where 



/sh(^)= / ds 



cos^(Z^^2^T/£) 
s^^/ s^ — 1 



(23) 



(24) 



(25) 



Below we will make use of the large Z limit of /sh(^), 
which is 



lim /sh(^) = t;- 



B. dz branches 



(26) 



For the to = + and — branches the normalized vibra- 
tional eigenmodes arei 



K r 



47rcA 1 



=F a 



^ — iaKz 



C± e'"-^") 



+ ^ 



iP^t 



-ifJKz 



+ C±e 



i0Kz\ 



eK 



I 4/ —iaKz I /- iaKz\ 

± a 2 ( e + C± e 1 



+ ir^(e-''"'^"-C±e*''^") 



e. >e 



iK- 



where 

a{c) = a/(c/wi)2 - 1. 

Here C± is a phase factor 

[(/32-l)±2zV^2 



C±(c) 



(/32-l)2+4a/3 ' 



(27) 



(28) 



with |C±(c)| — 1- We also note that 



^ei = ^(l-S...). 



(30) 



The m — ± branches therefore contribute an amount 



g'l+\r,^)+g'l^\r,Lu) 



(29) to the local DOS, where 
I 



(31) 



Here 



and 






(l-5..)/±(f)+2<5,.J±( 



UJZ \ 
V\ } 



I±{Z) 



'ds_ 



1 ' 



^— iZ\/i^^^/s 



T *(f^-l)^ 



-»Z^(s/,.)2-l/s _|_ 



(^ - 2)2 +4(,s2- 1)^(1^-2)^ 
[g;~2±2z(.^~l)i(g^-2)i]^ 



iZVs^-l/s 



jZ^(s/i.)2-l/s 






(.^-1)^ 



^-iZ\/j23l/s 



il2 



[^-2±2z(s2- 1)1(1, -2)3] 



^iZ\/iI^^/s 



± *^-l 



-iZ 



(£^_2)2 + 4(s2-l)^(^-2)i 

yilAO^A. (^ [^-2±2^(.^-l)n^-2)^]^ Y^^-^^/7IA^/ 
(^-2)2 + 4(.2-i)i(i^-2)i 

I 



(32) 



(33) 



(34) 



where 






(35) 



is the ratio of transverse and longitudinal bulk sound 
velocities. It can be shown that 



lim /± {Z)= , 



^ I ^ "^ (^ 2\3/2 ,„„^ 



and 



£m^J.(Z)^| + ^-|,(2 + .^)VT3^. (37) 



C. branch 

Next we consider the m — Q branch, for whichi 



fo 



K 



27r/3cA\ 



iCe-'^-^^ + i/3e-^'^^" 



where 



and 



i(3Ae 



i(3Kz 



Bk 



7C e 



-iKz 



le 



-ijSKz 



iAe 



ipKz 



AKt 



7(c) = ^1- (0/^1)2, 

(/32 - 1)2 - 4i(3-f 



Aic)^ 



(/?2-l)2 + 4i^7' 



C(c) = 4^(^' - ^) 



{13^ - l)2+4i/?7' 



(38) 

(39) 

(40) 

(41) 



Note that \A{c)\ = 1. 

The m = branch therefore contributes to Eq. p(i|) 
the quantity 



JO) 



g^'ir.u;) 



2^,0 



i>7T^v: 



(l-'5..)/o(^)+2(5..Jo(f) 



(42) 



where 



Io{Z) ^ 



Ijv 



ds 



s^V^ 



4(s2-2)(s2-l)i 



(s2_ 2)2 +4i(s2- 1)^(1 -jy2s2)^ 



^-ZVl-i'^sVs 



IPe 



— iZ^/s^ — 1/s 



+ {s^~iy 



(s2 - 2)2 + 4i(s2 - 1)^(1 _ ,y252)i 



(43) 



and 



MZ) 



Ijv 



ds 



1 s 



t\/i2~~T 



(1 



l.2s2)2 



4(s2-2)(s2_l)5 



(s2 _ 2)2 + 4i(s2 - 1)5 (1 - 1/2^2)^ 



^-ZVl-i^^s^/s 



le 



.zv^^/. , y (.^-2)2-4z(,s2-i)ni-^v)n ^.zv- 



(.S2 - 2)2 + 4i(s2 - l)i(l - iy2s2)t 



-lA 



(44) 



Furthermore, we note the asymptotic results 



hm /o(Z)--(l-z/2) 

zj — >oo o 



2n3/2 



and 



lim JoiZ) - - (2 + 1/2) Vl-:/2_ 



Z^oo 



(45) 



(46) 



where ^ is the root between and 1 of the polynomial 



f - 8^^ + 8(3 - 2i/2)^2 _ -^g^^ _ j^2) ^ Q^ (5^) 

with 1/ the velocity ratio defined above in Eq. H35|) . In 
terms of ^ we can write 



<y5 



VT~~^ and ?7 = v^l-C^- (52) 



D. Rayleigh branch 

The normalized vibrational eigenmodes for the 
Rayleigh branch m — H arei 



/ca\ 



le 



-(pKz _ ■/ 2y?) \ -(7K2 



1 1+772; 



(/je 



l+r;: 



Bk 



{j^)e-^'<^y^y^-, (47) 



where 



The Rayleigh branch contributes an amount 



where 



?lf^(r:w) 



AttJCcI 



(l-^-)^R(f) 



25..Jr(^) 



(53) 



if = v/1 - (cr/vi)^ and 77 = Vl - (cR/ft)2. (48) 
Furthermore, 

((/?-?7)(V?-77 + 2(/3772) 



/c 



2Lpvp 



(49) 



In Eq. (|48|l . cr is the velocity of the Rayleigh (surface) 
waves, given by 



and 



Ik{Z) 



-^z 



i^h-^f 



MZ) = ^-'[e 



= ,n2L-v^_ (_L-\^-vZ^ 



a+J72 



(54) 



(55) 



CR = ^Ut, 



(50) Both Ik{Z) and /r(2') vanish in the large Z limit. 



E. Local DOS 



By symmetry, the local DOS tensor is diagonal, with the in-plane components gxx — Qyy generally different than 
the perpendicular component g^^. Combining the results derived above, we obtain our principal result. 



gij{r,uj) = 5i 



2(2 + 1/3) 



Gtt^uj^ 



(l-<^..)5||(f)+'^..54f) 



where 



and 



5„(^)-^/sh(^/.) + ^^[/.(^)+/-(^)]+^(^/o(^/.) 



Stt 



2(2 + iy^)ICC- 



■iK{eie.i 



9±id) 



3v^ 



4(2 



[,/+(0)+J-(0)] 



2(2 + z/3) 



Jo{0/i 



Stt 



(2 + i/3)/c^: 



■jR(e/e^). 



(56) 



(57) 



(58) 



The prefactor u}'^{2 + i/^)/6tt'^v^ in Eq. H5()|) is the bulk Debye DOS. In the large 9 limit, corresponding to high 
frequencies u>, large perpendicular distances z, or both, the asymptotic results given above lead to 



lim gn{e) = lim gj_{e) == 1, 

9—>00 — >CX3 



(59) 



0.5, 



T 1 1 r 



g = a 




and mass density 
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FIG. 1: Dimensionless local DOS functions g\\ (solid curve) 
and g± (dashed curve), defined in Eqs. I|57^ and 1)58^ . for 
an isotropic elastic continuum model of a Si substrate. Here 
8 = ijz/vi. The DOS functions are essentially equal to unity 
for e > 10. 



in which case the local DOS of Eq. H56|l reduces to the 
Debye value given in Eq. H15|l . 



V. APPLICATION TO SILICON 

To apply these results, we shall approximate Si as an 
isotropic elastic continuum with longitudinal and trans- 
verse sound velocities 



vi = 8.47xl0^cms"\ 
vt = 5.86xl05cms-\ 



(60) 
(61) 



p = 2.33 g cm 



(62) 



The velocities above follow from using the measured elas- 
tic constants cn = 1.67 x 10^^ dyncm"^ and C44 = 
8.01x10^^ dyncm^^ at 70K reported in Ref. [3. A more 
precise method of determining these velocities would be 
to use the actual anisotropic elastic tensor to construct 
a best fit to an isotropic one.^ However this only makes 
small changes in the final results. 

The Rayleigh sound speed is determined from Eq. 151|l . 
The velocity ratio defined in Eq. H35|l is 1/ = 0.693, leading 
to ^ = 0.882, and hence 



cr = 5.17x10^ cm s 



(for Si) 



(63) 



In addition, ip = 0.792, 77 = 0.472, and /C = 0.610. 

In Fig. n] we present plots of the dimensionless DOS 
functions gu (0) and g±{0) for Si. Here = ujz/vi, where lo 
is the angular frequency, z is the perpendicular distance 
from the surface, or depth, and vi is the bulk longitudinal 
sound velocity. The dimensionless DOS has an interest- 
ing oscillatory structure when is small, and approaches 
unity in the large 6 limit. When = 0, 



and 



.911(0) = 2.60 



.g_L(0) = 3.28. 



(64) 



(65) 



On the substrate surface, z = 0, the local DOS is con- 
siderably modified from the bulk Debye value, as a con- 
sequence of Eqs. (|64|) and H65|) . In particular, g^xi^,^) 
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FIG. 2: Frequency dependence of the local DOS gxx{z,uj) — 
gyy{z,Lj) (solid curve) and gzziZyUj) (dashed curve) for a fixed 
z > 0, compared to the purely quadratic frequency depen- 
dence (dotted curve) of the Debye model. In all cases the 
DOS is given in units of (2-|-z/'^)/67r'^wjf . The frequency is plot- 
ted in units of vi/z, where vi is the bulk longitudinal sound 
velocity. Elastic parameters corresponding to Si are assumed. 



and gyy(0,uj) are quadratic in frequency with a constant 
of proportionality gy (0) times that of the Debye model. 
Similarly, gzz(0,uj) is quadratic in to, with a constant of 
proportionality a factor of g± (0) larger than in the Debye 
case. 

The oscillatory structure in g|| {6) and g±{d) also modi- 
fies the frequency dependence of the local DOS inside the 
substrate. This is illustrated in Fig.|21 where gij(z,uj) is 
plotted in units of {2 + i'^)/6tt'^v^. In these units, the 
Debye DOS is simply w^, shown as the dotted line in the 
figure. 



VI. DISCUSSION 



We have calculated the local, eigenfunction-weighted 
acoustic phonon DOS tensor for a model substrate con- 
sisting of a semi-infinite isotropic elastic continuum with 
a stress-free surface, and applied the results to Si. The 
local DOS on the surface is quadratic in lu, with a pro- 
portionality constant that is enhanced compared to the 
three-dimensional bulk Debye model. The enhancement 
factors for Si are given in Eqs. l|M)l and (|^ . Inside the 
substrate the frequency dependence is non-quadratic, as 
illustrated in Fig. |21 At depths z large compared with 
vi/u, where vi is the bulk longitudinal sound velocity, or 
at frequencies uj large compared with vi/z, the local DOS 
approaches the isotropic Debye value. These regimes cor- 
responds to the large 9 limits of gy {9) and g±{9), defined 
in Eqs. ^ and ifSS)) . 

As expected, the presence of the free surface signifi- 
cantly changes the vibrational properties of the system 
when the condition z < A is met, where A is a charac- 
teristic bulk phonon wavelength with frequency w. At a 
fixed frequency, this condition reduces to a requirement 
on distances from the surface, whereas at fixed depths it 
becomes a condition on frequency. 
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